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Abstract 

-^ We present a fast algorithm to calculate Coulomb/exchange integrals of 

prolate spheroidal electronic orbitals, which are the exact solutions of the 

single-electron, two-center Schrodinger equation for diatomic molecules. 

Our approach employs Neumann's expansion of the Coulomb repulsion 

i i 1/ | ec — j/l , solves the resulting integrals symbolically in closed form and 

_Ch subsequently performs a numeric Taylor expansion for efficiency. Thanks 

M^ to the general form of the integrals, the obtained coefficients are indepen- 

4— > dent of the particular wavefunctions and can thus be reused later. 

Key features of our algorithm include complete avoidance of numeric 
integration, drafting of the individual steps as fast matrix operations and 
high accuracy due to the exponential convergence of the expansions. 

Application to the diatomic molecules O2 and CO exemplifies the de- 
veloped methods, which can be relevant for a quantitative understanding 
of chemical bonds in general. 



> 

}£\ 1 Introduction 

The two-center electronic Schrodinger equation is a natural starting point to 
study diatomic molecules or chemical bonds. It is well known that it separates in 
~^ prolate spheroidal coordinates. Thus, the corresponding single-electron orbitals 

I can be calculated efficiently. For several electrons, however, the tedious inter- 

electron Coulomb repulsion integrals have impeded a widespread use of these 
orbitals so far. To alleviate these difficulties, we present an efficient algorithmic 
framework in this paper. 



In the computational chemistry literature, LCAO (linear combination of 
atomic orbitals) is the most common approach to construct electronic wave- 
functions for molecules. It dates back to the early days of quantum mechan- 
ics [21]. In the seminal paper [5], Boys proposed Gaussian-type atomic orbitals 
since the necessary integrals can be explicitly evaluated. Hence they are widely 
used in modern computational chemistry software packages. Nevertheless, only 
the exact single-electron spheroidal orbitals are - by definition - precise for any 
distance of the atomic nuclei. This fact is an important advantage for studying 
diatomic molecules and chemical bonds. 
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An interesting alternative approach to diatomic molecules is the Holstein- 
Herring method [16, 14, 33, 30] for calculating exchange energies of H^-like 
molecular ions. This method has recently been extended to two-active-electron 
systems [28]. However, it is not suitable for an arbitrary number of valence 
electrons. 

Another common approach first proposed by Hylleraas [18] for the helium 
atom includes the inter-electron distance r^ as independent variable into the 
electronic wavefunction. Thus, the pairwise electronic Coulomb cusp is han- 
dled explicitly, which reduces the number of required wavefunctions. James 
and Coolidge [19] have applied this method to the H 2 molecule using spheroidal 
coordinates, which still serves as starting point for modern benchmark calcu- 
lations. Ref. [7] contains an extension to the HeJ and He2 molecule, and a 
modern review can be found in [20]. 

Ref. [9] is part of a series which provides an extensive discussion of Gaussian 
basis sets for molecular calculations, and specifically computes the total energy 
and dissociation energy of O2. 

Ref. [23] employs Kohn-Sham density functional theory for diatomic molecules 
in (discretized) spheroidal coordinates. In particular, the authors apply their 
method to calculate the ground state energy of carbon monoxide CO. 

The basic setup of prolate spheroidal orbitals employed in the current paper 
has been developed in Ref. [2, 3] and applied to molecules with up 4 electrons. 
Our contribution is a reformulation into an efficient computational framework 1 , 
which allows for an extension to many more electrons. For example, the oxygen 
dimer O2 contains 16 electrons. 

Outline Section 2 provides the details of the single-electron Schrodinger equa- 
tion for atomic dimers in prolate spheroidal coordinates. Our presentation is 
based on the series [2, 3], and additionally includes a "best match" mapping 
to the common LCAO molecular orbitals. Section 3 contains the main ab- 
stract mathematical contribution of this paper: we prove a recurrence relation 
to efficiently multiply function expansions in terms of associated Laguerre poly- 
nomials, and solve several integrals symbolically in closed form. These results 
(combined with Neumann's expansion of 1/ \x — y\ into Legendre polynomials) 
are the basis of our algorithm. It is described in detail in section 4, including 
cost analysis and error estimation. Section 5 contains the application of the 
algorithm to the O2 molecule, which is particularly interesting among atomic 
dimers due to its paramagnetism. 

2 Single-Electron Schrodinger Equation for 
Atomic Dimers 

This section introduces the single-electron quantum mechanical framework, which 
serves as starting point for the many-electron calculations in section 4. We ba- 
sically follow the discussion in Ref. [2, 3]. 



lr The complete source code of our implementation is available online at [26] (in the 
mathematica/diatomic subfoldcr) 



Separation in prolate spheroidal coordinates The single-electron, two- 
center Schrodinger equation for a ITJ -like molecular ion in atomic units (Born- 
Oppcnhcimcr approximation) reads 



2 r a r b 



(2.1) 



Here, r a and r& denote the distances to the fixed nuclei at (0,0, =pR/2), respec- 
tively, and Z a ,Zb G N>o the nuclear charges (see figure 1). The distance R 
between the nuclei is also called bond length in the chemistry literature. We 
have omitted the repulsive interaction of the nuclei ( a R b ) for now to focus on 
the electronic energy, but will include it into the total energy later. In what 
follows, we set Z := ^±^ and Aq := (Z a - Z b ) R (w.l.o.g. Aq > 0). The 
homonuclear case corresponds precisely to Z a = Z b = Z and Aq = 0. 



z A 
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Figure 1: Spatial arrangement of a single electron bound to two atomic nuclei 



It is well known that equation (2.1) is separable in prolate spheroidal coor- 
dinates (£,r],<p) defined by 

£ := (r„ + n) /R, e > 1 
V-{r a -r b )/R, ??e[-l,l] 



and the Ansatz 



lH6»?.¥>)=A(OS(T7) 



s/2i 



(2.2) 



to G It is the eigenvalue of the angular momentum operator L z — —id^ , which 
commutes with the Hamiltonian on the left hand side of (2.1) due to the az- 
imuthal symmetry about the internuclear axis. In the following, we set /i := \m\ 
to shorten notation. 

Plugging (2.2) into (2.1) leads to coupled ODEs for the radial part A(£) and 
angular part S(rf). The latter reads 



d_ 



((l-^ 2 )^)+(^ J 4)-A,r ? + (i P ) 2 (l-r ? 2 )- r 



A?(ip,Ag) 



S(V) = o, 

(2.3) 



where the "energy parameter" p £ R>o is defined via the energy E, 

E=:-2(p/R) 2 , (2.4) 

and A is an eigenvalue of the operator Q (defined in [2]). For the purpose of 
this paper, we simply regard A as separation constant. In the homonuclear case 
Aq = 0, equation (2.3) is well know as the angular spheroidal wave equation [32, 
24, 11, 10] when we identify (p 2 — A) as spheroidal eigenvalue X^(ip). Successive 
£ — /i, [i + 1, . . . label the discrete set of eigenvalues for which (2.3) has a 
normalizable solution. 

Since E is finite, p —t in the united atom limit R — > 0, and (2.3) reduces to 
Legendre's differential equation. Then limjj_>.o Q = —L (angular momentum 
operator) with eigenvalue A = — A^(0) = —£(£ + 1). However, except for this 
special case, I is no valid quantum number since L does not commute with the 
Hamiltonian in general. 

The homonuclear solution S(ri) = S^(ip,r]) is already built into Mathemat- 
ica 2 and could thus be plugged into (numeric) integrals. Nevertheless, in order 
to use some properties of Legendre polynomials later and cover the heteronuclear 
case also, we employ the series expansion 

S?(ip,Aq, V ) = f^4 tk {p,Aq) ( 2fc + 1 j* ~ ffl * Kin). (2-5) 

Plugged into (2.3) results in a three-term recurrence relation for the coeffi- 
cients c^ u(jp) = c^ k (p,0) (homonuclear) and a five-term recurrence relation for 
Cp k (p,Aq) (heteronuclear) [24, 2]. Namely, in the homonuclear case, only inte- 
gers k with the same parity as £ contribute to the sum due to symmetry. After 
truncating this expansion (which is justified due to the exponential decay of the 
coefficients), it may be rewritten as eigenvalue equation (see also [15]) 



F»(p 1 Aq)c = \c, c=[c% tk (p,Aq)) k , A = A^(ip,Ag) (2.6) 

with a symmetric matrix F^f^p, Aq). This matrix is tridiagonal in the homonu- 
clear case (after proper relabeling) and pentadiagonal in the heteronuclear case. 
Note that fast eigenvalue solvers exist particularly for tridiagonal matrices. We 
adopt the normalization scheme used by [24] and Mathematica, namely 

j\(ip,A q , V fd V = JT \cl k (p,Aq)\ 2 I 2 / +1 [^j . (2.7) 

The energy parameter p couples (2.3) to the radial equation 



9 f (e-i)^)-(P 2 -A)+2ZRt+(ip) 2 (e ^ r 



«9£ V 9U ' i — -s—s £ 2 -l 

A£(ip,Ag) 



A(0 = 0. 



(2.8) 



2 Specifically, the implementation [10] has been integrated into Mathematica as 
SpheroidalPS[n,m,7,z] and SpheroidalQS[n,m,7,z] for the angular spheroidal function of 
the first and second kind, respectively. 



This is the radial spheroidal differential equation except for the 2 ZR £ term, 
and formally resembles (2.3) apart from £ > 1 versus \rj\ < 1. We determine p 
numerically as follows. 

First, define Hylleraas functions via associated Laguerre polynomials as 



H£(x) 



— T C/2 P -V2, 



./- -<■ -' 2 y/kl/(k + n)\ L%(x), k, fie N . 



This choice precisely incorporates the orthogonality relation for Laguerre poly- 
nomials, such that 

/>oo 

(2.9) 



(2.10) 



H£,(x)H£(x)dx = 6 k k>. 

Given a sequence d = (dk)k>o> we se t 

oo 

fc=0 



(Note that k starts at instead of \i as in (2.5).) Employing such an expansion 
for the radial wavefunction, 



A(e) = F£(2p(£-l)) 



(2.11) 



results in a three-term recursion formula [2] for the to-be determined coefficients 
dfc. They will turn out to decay exponentially, as illustrated in figure 2. Hence 
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Figure 2: Laguerre expansion coefficients of the (£,m) = (0,0) groundstate 
radial wavefunction (see equations (2.10) and (2.11)). The exponential decay of 
the coefficients justifies the truncation of the expansion. 

we can truncate the expansion and rewrite the recurrence relation as matrix 
equation [2] 



(B»(p)R"(j>, A) + Pf i 2 l) di 0, A = \%(ip, Aq). 



(2.12) 



Both R^ and B^ are symmetric tridiagonal matrices, and / denotes the identity 
matrix. The left hand side is singular for a discrete set of values p only. This 



condition finally determines p and the energy E. Ref. [2] employs a Newton 
iteration to obtain both p and A simultaneously, such that the matrices in (2.12) 
and F^ — X I in (2.6) have zero determinants. An improved version uses the 
so-called Killingbeck method [13, 29]. In our case, we apply a numerical root 
search algorithm over p such that an eigenvalue of the matrix in (2.12) becomes 
zero. 

Considering the starting point of the numerical iteration, [2] usespo = ZR/n, 
which becomes exact in the unified atom limit R — > and is thus valid for 
small ZR. Here, n labels successive eigenvalues as in the unified atom limit. 
Alternatively, we have identified p = ZR/(2n) as reliable candidate for large 
values of ZR, which stems from the dissociation limit R — > oo (hydrogen-like 
atom plus isolated nucleus). 

Figure 3 shows the lowest few homonuclear energy levels in dependence of 
ZR, both with and without the (rescaled) nuclear repulsion term \j(ZR). In 
analogy to the molecular term symbol, we employ the notation 

nt 2s+1 m g/u (2.13) 

to label states. In common notation, £ = 0, 1, 2, 3, . . . is designated by s, p, d, f, . . .. 
respectively, and m = 0, ±1, ±2, ... by a, ±7r, ±<5, . . .. For fixed (£, m), the "prin- 
cipal value" n — 1,2, .. . enumerates successive energy levels. In the homonu- 
clear case, the angular spheroidal wave function determines the parity (— 1) £ (re- 
flection about the origin, x — > —x). It is written as gerade (even) or ttngerade 
(odd). We omit the spin variable s for now, which will become important for 
the many-electron calculations in section 4. 

Having the exact solution of the two-center Schrodinger equation available 
calls for a comparison with the popular LCAO approach (linear combination 
of atomic orbitals). Figure 4 tries to match the corresponding wavefunctions, 
taking parity and ordering of energy levels into account. However, note that 
the suggestive ordering has to be interpreted with caution since it depends on 
the nuclear distance R. For example, according to figure 3, 

-E'lso-g < -^lpcr u j^lp(±7r) u j^2scr g < -E^po^ , ■ • ■ 

for small nuclear distances R. This is different from the arrangement in figure 4. 

Normalization In what follows, we derive a formula for the required normal- 
ization factor of the wavefunction. The volume element in prolate spheroidal 
coordinates equals dV^ = (-R/2) 3 (£ 2 — rj 2 ) d£ d?? dip. Thus 

/OO />1 
J ^A(0 2 S( V ) 2 (e-V 2 )d V dt 

The inner integral without the factor rj 2 is already solved in (2.7). To include 
rj 2 , we use the identity 

* • W - ^n-M + k ~^n + M- (2-14) 

Thus, after taking into account the normalization factors in the expansion (2.5), 
we obtain 

^S( V ) 2 V 2 d v =\\xi: cg c\\ 2 , c=(c£ fc (p)) fc , 
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Figure 3: Single-electron energy levels (equation (2.4)) of a Hj~-like homonuclear 
dimer with respect to ZR (nuclear charge x nuclear distance), in atomic units. 
The bottom plot additionally includes the rescaled nuclear-nuclear repulsion 
term 1/(ZR). The unified-atom limit R — ¥ corresponds to a hydrogen-like 
atomic ion with one electron, nuclear charge 2Z and energy levels ~2Z 2 /n 2 , in 
agreement with the curves of the top subfigure. In the dissociation limit R — > oo, 
the dimer splits into a single hydrogen-like atom/ion and an isolated nucleus 
(H + p for Z = 1). Thus, the electronic energy levels converge to — \Z 2 jr? . 



with the tridiagonal, symmetric matrix X? given by 



xt. u, = o, X M 



Lce.fcA; 



Lcg,fc,fc+1 



(fe + l-^)(fc + l + M) 

(2fc + l)(2fc + 3) 



1/2 



k = n, n + 1, 



We proceed analogously for the radial part. After a change of variables 
x := 2p(£ — 1) and due to the orthogonality (2.9), we obtain 



/ A(0 a d* = -||d|| a , 



2p z <-> lftTu 




2p 2 •<-> ldo-g 



2px(^2p y 2p z 



2s -@- 




2s -H- 2pcr u 



(T g 2s <-> 2sc g 



-@-2s 



i fl -e. 



Is <-> lpfT u 




IS f")- lSfTp; 



-e-i s 



Figure 4: Putative best match of the exact H^~-like electronic wavefunctions 
(labeled n£m g / u in boldface blue) with the LCAO-MOs (molecular orbitals built 
from linear combinations of atomic orbitals) widely used in the literature (see 
e.g. [1]). In particular, the parity (reflection about the origin, x —> —x) agrees 
in each instance. Orbitals are schematically drawn in red, and antibonding MOs 
are marked by a star (*). 



where d contains the expansion coefficients in (2.11). To incorporate the factor 
£ 2 , we employ the following well-known relation for Laguerre polynomials: 

x ■ L*(x) = -(k + l)L% +1 (x) + (2k + ft + l)L^(x) -(k + i*)L*_i(.x). (2.15) 

Thus, multiplying an expansion (2.10) by x yields 

x ■ H%(x) = H%,{x), d'—X^d (2.16) 



with the tridiagonal, symmetric matrix X£ defined by 



*La g)fc fe = 2fc + M + l, 



^^ fe . fe+1 = -((fc + l)(fc + M +l)) 1/2 



Lag,fc,fc+ 

Plugging (2.16) into the following integral yields 



f a<ov*-£ 



/ + (2p)- 1 X£L )d 



fc = 0,l, 



Assembling the radial and angular contributions finally results in 



Ml? 



(R/2f 
2p 



d+(2p)-'X£ a d 



X 



Leg 1 



That is, we obtain the correct normalization factor directly from the expansion 
coefficients c and d. 

Dissociation limit R — > oo From a physical point of view, separating the 
nuclei from each other should yield a hydrogen-like atom/ion plus an isolated 
nucleus. However, in the homonuclear case, the symmetry properties of the elec- 
tronic wavefunctions (ip(x) = (— l) e ip(— x) due to parity) imply that the elec- 
tronic charge is equally distributed to both nuclei. This seeming contradiction 
can be reconciled by constructing superpositions of even and odd wavefunctions 
to obtain the well-known hydrogen-like wavefunctions, localized at either one 
or the other nucleus. (Note that conversely, the LCAO approach uses linear 
combinations of atomic orbitals as molecular wavefunctions.) 

From the above arguments, we expect the energy levels to converge to 
— 2^ 2 1™ 2 m the limit R — > oo, as indicated in figure 3. Along with it comes 
a heuristic understanding of the convergence rate 3 . Each "half" electron local- 
ized at a nucleus experiences an additional attraction from the respective other 
nucleus. This adds up to the net attraction energy 



0.5 x Z 
R 



0.5 x Z _ Z 
R ~ ~R' 



(2.17) 



Subtracting this correction term (which of course vanishes as R — > oo) from the 
energy E leads to exponential (instead of algebraic) convergence, as shown in 
figure 5. Namely, the electronic charge distributions decay exponentially with 
distance from the nuclei, implying a likewise decay of the error. 



|£u,r/Z z +l/(ZR)+l/2| 

io- 4 - 




ZR[a.u.J 



Figure 5: Exponential convergence of the lscr energy level E/Z 2 plus the ^ 
correction term (2.17) to the hydrogen groundstate energy — |, as R — > oo. 



3 1 am grateful to Gero Friesecke for helpful discussion regarding this point. 



3 Properties of Laguerre Expansions 

This technical section is based on function expansions in terms of associated 
Laguerre polynomials (see equation (2.10)). We develop a computational frame- 
work for multiplying these expansions, and derive analytic solutions of integrals 
appearing in section 4. 

3.1 Products of Laguerre Expansions 

We want to solve the following task: Given integers morris € Z and exponen- 
tially decaying sequences (di^), (d,2 t k)> calculate the sequence (df~) satisfying 

H^\x)-H^{x)=H^- m2 \x). (3.1) 

For conciseness of notation, let [ii :— |toj|, i = 1,2 and /U3 := \m\ — m®], and 
assume without loss of generality that [i\ > /X2- Depending on the signs of 
rtii and 7712, we have /13 = \X\ ± \Xi- The orthogonality relation of Laguerre 
polynomials leads to 

d fe = <d 2 |n£di), k = 0,1,2,... (3.2) 

with the symmetric matrix LT^ = (afj k )ij given by 

Hg(x)Hg(x)Hg(x)dx 







1 



1/2 ( „ (3.3) 

/ 6f (3/2) if M 3 = Mi- 



tt »fc! \ „ J °i l d / 2 J ir A*3 = Ml - M2 



fi 1 (4+Mfe)!y \^(3/2) if M3 = /ii+M2 

In the above expression, 

/*OC 

&f(z):=/ /'^^^(^(xle-" da;, (3.4) 

JO 

^0*):=/ a^ +w ££(aOl£ (*)££(*) 6"" dx (3.5) 



/o 
defined for z G R>o- Using the recurrence relation 

x»L?(x) = {i + »)x»- l Lr l (x) - (i + l)*"" 1 !^*), 

the integrals (3.4) and (3.5) can be reduced to the following proposition, which 
is a generalization of [22] . 

Proposition 1. Given fixed integers /x £ Nq, the coefficients 



cf(z):= L%(x)L%(x)L%(x)e-"iix, z£l>„ (3.6) 

Jo 

defined for i £ N§ obey the recurrence relation 
c?(z) = - (1/z - 1) K_ lii2 , 3 (z) + <, fa _ lii3 W + cr ili2 ,3-iW) 

+ (V* - 1) Kfa-Ma-lC*) + <-l, 42 ,3-l(^) + <-l,i,-M,(*)) 

-P/^lX-^-r.-iW (3.7) 

10 



with the convention that c?(z) = if any i\, i 2 , i 3 < and r . j = doi for integer 
i > 0. 

Thus, c^(z) can iteratively be calculated and stored for later usage. Note that 
the coefficients c?(z) are symmetric in i\,i 2 ,i 3 . The case z = 1 and /x = 
is handled in [22] (with a sign typo in his equation (7)). From the particular 
form of the binomial coefficients in (3.7) it follows that the recurrence relation 
is homogeneous precisely if any fi k = 0. 

Proof. A derivation of (3.7) proceeds along the same lines as in [22], involving 
generating functions of Laguerre polynomials. More specifically, using 

°° e xt/(i+t) 

4=0 V ' 

the following formal series in ti,t 2 ,t 3 fulfills 

G»(t t ,t 2 ,t 3 ,z):= £ c^z)f[(-t k )^=z- ^=f + tk) ~ n (3.8) 
H,i^=o fe=i l-p(ti,h,h,z) 

with 

p(h,t 2 ,t 3 , z) ■- (i/z - 1) (tj + t 2 + 1 3 ) 

+ (2/z - 1) (ht 2 + ht 3 + tats) + (3/2 - l)tit 2 t 3 . 

Applying the identity 1/(1 — a;) = 1 + x/(l — x) for x — p(t\,t 2 ,t 3 , z) to the 
right hand side of (3.8) leads to 



1 3 

G»{t 1 ,t 2 ,t 3 ,z)= P {t 1 ,t 2 ,t 3 ,z)G>*(t 1 ,t 2 ,t 3 ,z) + -~[[(l + t k )-^. 



fc=i 

Now comparing coefficients of t^ t l 2 2 t 3 3 gives equation (3.7). □ 

Numeric experimentation suggests that c®(z) is bounded asymptotically 
(|i| — > 00) if and only if z > 3/2. As illustration, figure 6 shows the central 
coefficient c°,-j(3/2), which alternates its sign depending on the parity of i. 
As heuristic explanation of the asymptotic behavior, we focus on the central 
coefficient c° t ^z) and set 

c 3l (z) := cf )M («), Csi-i(^) := c°_ Mji (z), £3,-2(2) := c°_ m _ m (;j). 

Plugged into (3.7) and letting k := 3i gives 

c fe (z) = -3 (1/z - 1) 4-1(2) + 3 (2/z - 1) c k - 2 (z) - (3/z - 1) c k _ 3 {z). 

This equation is only correct if A: is a multiple of 3. Nevertheless, interpreted as 
difference equation yields the companion matrix 

/ 1 

z= 1 

\-(3/z-l) 3(2/z-l) -3(1/^-1), 



11 



4(3/2) 



0.05 



60 .80 J00 



Figure 6: Asymptotic behavior of the central coefficient c° ii (3/2) defined 
in (3.6), which oscillates between positive and negative values. 

with eigenvalues {1, 1, — (3/z— f )}. Thus, the spectral radius p(Z) < 1 precisely 
if z > 3/2. 

We add the following observation: the homogeneous recurrence relation may 
be interpreted as a differential equation in 3 dimensions by treating the indices i 
as continuous variables, cf(z) = f^{i), and taking the continuity limit. Namely, 
without the inhomogeneous contribution, (3.7) becomes 

= p[ - (1/z- i) (Jit d - (M,o)) + #(» - (0,h,0)) + f?(i - (0,0, h))) 
+ (2/z - 1) (f?(i - (0, h, h)) + f?(i - (h, 0, h)) + /**(t - (h, h, 0))) 
-(3/z-l)f?(i-(h,h,h))-f?(ij 

= -\ (d t2l3 + d ili3 + d lll2 ) #(») + h (J- z - l) 9 ul2l3 #(t) + O (h 2 ) . 

Here we have already used 

(d i2h +d ili3 +d ili2 )f?(i)=0 
to simplify the 0(h) term, which disappears precisely for z = 3/2. 

3.2 Argument Rescaling 

Given any fixed y € R>o, we try to re-express Laguerre expansions (2.10) eval- 
uated at the scaled coordinates yx as expansions evaluated at x. First note the 
following well-known identity for fc,/i£ No, 

L£(yx)=y k J2( 1 /y- 1 ) k - i ( k t'*) L i(x) for y e R >0 . (3.9) 
Similarly, a direct calculation shows that for all y ^ 0, 

^(i- 2 /r'xr(^)=/'^(i/ y -i)^(J+' i + 1 )Lr(x). ( 3 .io) 

i=0 i=0 ^ ' ' 

12 



Due to (3.9), for any exponentially decaying sequence d := (dk)k>o it holds that 
H${yx) = y^e-^-W 2 H£(x), d! := S»d (3.11) 

with the upper triangular matrix S^ = (sf fe (y)) defined by 



k\ (t + ^y.\ 1/2 k ^fe-iA + M 



«£(*):= (££S) y*(V» "I)*"'! ::?) ^ ,<fc 



and s^ k {y) = otherwise. 

This result can be combined with the operation (3.1), as follows. Assume 
we are given /x € Nq with /13 = jj,\ ± fj,^, as well as 21,22 € K>o and two 
exponentially decaying sequences d\,d2. Set z := (zi + zi)j2 and use (3.2) to 
calculate the sequence d = (dk)k>o, 

d k := (zir /2 (4r /2 (sg d 2 I Il£ S% di) , z< := %/*. 

Then, combining (3.1) with (3.11) gives 

H£(z 1 x)H£(z 2 x) = H^(zx). (3.12) 

Note that the exponential functions on both sides match. Summarizing, we have 
obtained the product of two Laguerre expansions with rescaled arguments. 

As slight variation of (3.11), given y € M>o and an exponentially decaying 
sequence d, we try to find a sequence d such that 

H^{yx) = H%{x). (3.13) 

Due to the orthogonality relation (2.9), we have to compute the following inte- 
gral for integers i, k > 0. Using equations (3.9) and (2.9) leads to 

/ -. \ i+k / . 



x (~ 1 ) ^-tt 2^1 -*,-fc;l + w- 



with the Gaussian hypcrgeometric function 2 -Pi- 

3.3 Integral Identities 

We derive analytic solutions of integrals originating from Neumann's expansion 
of 1/ \x — y\ in terms of Legendre functions (see equation (4.2) below). 

Proposition 2. For any y,z G R>o and integers fc,/i > 0, it holds that 

L»(yx)e- x dx=(l-y) k -L£(yz)e- z 
fc-i / ,. . x x (3-14) 



?=o 



^(l-r i - i (l/if(F)e-+(^ 
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Proof. First note that for y = 1, equation (3.14) simplifies to 1 — c z for k = 
and 

r^(x)e-dx^(L^ 1 (z)-^(z))c-+^ + ^ 1 ) for k > 1. (3.15) 



This identity can be proven by taking derivatives on both sides. Then, combin- 
ing (3.9) and (3.10) (for k - 1) with (3.15) leads to (3.14). □ 

For any integers /i, k, k > 0, consider the nested integrals 

^ := 2./ qWe- £/2 J o LZ(x)e-*/ 2 dxdx. (3.16) 

They have a surprisingly simple form for /j, = 0, 1: 
Proposition 3. TTie integrals t - in (3.16) are equal to 



kk 

1 k = k 



'1-2 2-2 

ti k = {2 { -i)^ k< k (4) = 1 1 ill 

fc > fc 



and ifte integrals t - 

{{— l) k k <k and k even 

1 k > k and k even v-hi) = 

otherwise 



/ 1 -l l-i 
/ 1 o o o 



10 1-1 
10 10 



V 



Proof. These identities can be proven by applying proposition 2 to the inner 
integral and using the orthogonality property of the Laguerre polynomials. □ 

For the following paragraph, we state 
Definition 4. Given integers /j£No and 1 < i < k, set 



Expressed in terms of generalized hypergeometric functions, 

'k + (j\ ( 1 1 i-k 

J + fi) 3 2 {l + i 1 + i + fJ- 
Given z € K>0j 2/ £ K and integers fc, /i > 0, we set out to solve the integral 

L£(ya;)e- x log(l + -) da. (3.18) 

For that purpose, we decompose the logarithm into log(l + x/z) — log(;r/2:). 
Considering the first term, integration by parts and (3.14) give 

f°° 1 
log(l + x/z)I%(yx)e- x dx = / L£(yx)dx 

Jq x -\- z 

-y^il-y)"- 1 - 1 L^yx)e~'dx. 

l=0 Jo x + z 

The integrals on the right hand side are solved by the following proposition: 
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Proposition 5. Let z £ M>o and y £ M., then for all integers k, \i > ; 

^oo -i , k 

10 ^ - - i=1 

with the incomplete gamma function T. 

Proof. For k — 0, we obtain (via a computer algebra system) 



/ -—L£(y X )c-*dx = LZ(-yz)r(0,z)c*--J2KM (JL r- (3-19) 
J x + z z f— ' i! 



,00 1 

/ e _x da; = r(0, z) e z 

Jo x + z 



in agreement with the right hand side of (3.19). For k > 1, a change of variables 
yields 

-^- (yxrc^d.T = f° -L- (yz.fe-d, - (_ yz )^r(0,z)e z 
a; + z J x + 1 dz % 

for any integer i > 0. Thus, the integral (3.19) is a linear combination of the 
last term (2 = 0, . . . , k). The explicit formula on the right hand side of (3.19) 
follows from a (rather tedious) calculation, using d z T(0, z) = —cT 7 -jz. O 

Concerning the second logarithm log(x/z) in the above decomposition, first 
note that 

f°° x { 

fi(z) := / log(a:/z) — c x dx = H t - (7 + log(z)), 
Jo * ! 

where 7 is Euler's constant and Hi the i th Harmonic number. Namely, integra- 
tion by parts yields the recurrence relation 

.n(z) = \ + f l - 1 {z), i = i,2,..., 

and fo(z) = — (7 + log(z)) can be shown by a computer algebra system. Thus, 
for all integers k, [i > 0, 

J™ log(x/z)l£(x) e"* dx = J2 (*+£) ("im -[ k + I *) (7 + log(z)). 

Combining this equation with (3.9) yields the following generalization: 
Proposition 6. Let z € K>o and t/6R, i/ien for all integers k, fj, > 0, 



f°° log (x/z) Lfryx) e"* dx = £ f fc + M ) (-y)'F, 
Jo i=1 V* ■+■ M/ 



k + Afn Mi \r^{k\y i (l-y) k - i * 



k ) \ r-' \i/ i + /i 



i=i 



(l-g) fc +A«E / (7 + log(,)) 



Hence we have collected all ingredients for solving the integral (3.18) in 
closed form. 
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We can now assemble the above results to calculate the nested integral 

L£(£)e- £/2 arcoth(l + £/z) / L%(x)e~ x / 2 dxdx, zeR >0 (3.20) 



for integers k,k,/i <G No- Proposition 2 with y = 2 and a change of variables 
(2x — > x) gives the inner integral. Combined with the Laguerre product coeffi- 
cients cf (1) in (3.6), 

(3.20) = 2 q^ k [ Vipx) e~ x log (l + -) da; 

k+k „oo , „ \ 

-E<^/ ^Wc-log 1 + - dx 
with the integer (!) coefficients 

< fc =-(-i) fc +D- 1 ) fc - 1 - i (*!?)• 

i=j-k 

The two above integrals are precisely of the form (3.18), which completes the 
calculation of (3.20). 

As last task of this section, given k, /i € No and z G K>0j we try to compute 
the Laguerre expansion coefficients of 



V?z + 



In other words, due to the orthogonality property of Laguerre polynomials, we 
have to calculate the integrals 

( X . 1 H£(x)W{x)dx (3.21) 

7o V2z + x k 



the above integral can be reduced to a linear combination of 



for k = 0, 1, Employing the Laguerre product coefficients b^. (1) from (3.4), 



w k (z):= L k (x)e- X dx, z£l>«. (3.22) 

Jo \J2z + x 

We can calculate these integrals iteratively for k = 0, 1, 2, . . . via the following 

Proposition 7. 77ie functions Wk{z) obey the recurrence relation 

Wfc(z) = Wfe-i(z) + — — ■= — , fc = l,2,... (3.23) 

fc ax y z 

wit/i i/ie starting value 

wq(z) — v^erfc ( \Tlz\ e 2z , 
where erfc is £/ie complementary error function. 



1G 



Proof. The formula for Wq{z) can be derived via a computer algebra system. 
Concerning the recurrence relation, first integrate (3.22) by parts (x — > x/z) to 
obtain the alternative representation 

Wk(z) = y/z , L k (zx)e zx dx. 

J y/2 + x 

Applying the relation 

^-L k (zx)e-" = t±±(L h+1 (zx)-L k (zx))e-' x (3.24) 

dz z 

to this representation gives the recurrence formula (3.23). The relation (3.24) 
follows from combining 

d fc 

— L k {z) e~ z = - ^2 L i( z ) e ~ Z with 

i=0 

x ■ Li(x) = -(i + l)L l+1 (x) + (2i + l)Li(x) -iL^x). 

D 

4 Coulomb and Exchange Integrals of Prolate 
Spheroidal Orbitals 

The computation of Coulomb interactions is often the most demanding task 
concerning multi-electron quantum systems. In this section, we provide the 
details of an efficient algorithmic implementation, which employs analytically 
precomputed integrals (from section 3) and a subsequent Taylor expansion to 
speed up calculations, and avoids difficulties caused by an alternative numeric 
approach. For example, we observe that the absolute value of the nested inte- 
grals in equation (4.9) is typically much smaller than (the maximum over x) of 
the inner integral. This general effect could be explained by the orthogonality 
property of Laguerre polynomials. In any case, analytically solving the nested 
integrals as a whole circumvents the numeric difficulties caused by the blow-up 
of the inner integral. 

Given square-integrable spatial "orbitals" a,b,c,d € L 2 (R 3 ,C), we define 
the Coulomb integral (following standard notation) as 

1 



(ab\cd):= / a(xi)b(xi) -. r c(x 2 )d(x 2 ) dx ± x 2 , 

Jrb \Xl-X 2 \ 

where ~ is the complex conjugation. In our setting, we want to calculate the 
concrete realization 

(^ntm^n'i'm' l^nl-m^u'l'm') ( 41 ) 

for single-electron wavefunctions ip n em from section 2. The labels n£m, n'£'m' 
etc are the "quantum numbers" in the molecular term symbol (2.13). 

To evaluate these Coulomb integrals in prolate spheroidal coordinates, we 
pursue the same approach as [2] and employ Neumann's expansion 

1 4^^. _ 2r + l /(r-*/)! x:i 



|Irtrj i:D-.r^(^)™« (42) 

xP^O^Mcos^i-^)) 
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with eo = 1, t v = 2 for v > and £i < £2 (otherwise interchange £1 -H- £2)- P" 
and Q V T are the Legendre functions of the first and second kind, respectively. A 
derivation of (4.2) can be found in [27]. For the following, remember the volume 
element in prolate spheroidal coordinates, dV^ = (-R/2) 3 (£ 2 — -q 2 ) d£d7/dvJ. 
With (4.2) plugged into (4.1), the integrals over ipi and (p 2 result in 



1 



2-7T f 2lT 



cos(r/(<^ - tp 2 )) e -Kv>-™')vi-i<fr-*')v* di Pl dip 2 



(27r) 2 7 Jo 

I l/e„ if to — m! = — (ttt, — to') and r/ = |m — to'| 







otherwise 



(4.3) 



Thus, the right hand side of (4.2) effectively contains the sum over r only, 
starting from t = v. 

To approximate the infinite sum over r, we include all terms up to a thresh- 
old r max . This truncation is justified due to the exponential convergence, as 
illustrated in figure 7. 



0.001 



Up* lp(-ff) I lp(-;r) lpi) 




(lscr lscr | lstr lscr) 
(lscr lptr I lpcr lscr) 



Figure 7: Estimated relative truncation error of the sum over t — v, . . . , r max 
in the Neumann expansion (4.2), exemplified by the Coulomb integrals (4.1). 
Wavefunctions are taken from section 2 with R = 121 pm (experimental bond 
length of 16 C>2). The observed exponential convergence renders the expan- 
sion (4.2) particularly useful. We have calculated the error by comparison with 
the sum up to r = 11. Only every second r contributes due to the symmetry 
constraint (4.5) below, which explains the plateaus of the curves. 



4.1 Angular Coulomb Integral 

The expansion (4.2) admits a separation of the r\\ and r\i integrals. Both are 
of the same form, so it suffices to restrict the following presentation to the r\\ 
integral. Taking into account the volume clement, we have to calculate 



S?(ip, Aq, V )sg (ip', Aq, 77)P;(7 7 ) rf drj 



(4.4) 
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for j G {0,2}, where we have once again set fi := \m\ and pi := \m'\. Since 
S%(ip, 0,77) and P"(v) nave parity (— l) e+tl and (— 1) T+ " , respectively, it follows 
that in the homonuclear case, (4.4) is non-zero only if 

t + fi + l' + n' + T + v is even. (4.5) 

Plugging the expansion (2.5) into (4.4) results in a linear combination of inte- 
grals of the following form (see also [2, Appendix D]), which arc explicitly solved 
by Wigner 3j symbols: 





f{ti-IH)\\ 


f PtiWPgWP&Wtv 








'(_!)»+» (>' >' >° 


= 


fix l 2 V 

K o 0, 


)*. 


\Ml ~(J>2 ~M3 
VMl M2 -M3/ 



if ^3 = |/ii - p 2 \ 



if A*3 = A*i + /"2 

The equation is valid for non-negative integers p,\,p2,P'3, assuming (w.l.o.g.) 
Ml ^ M2- Given the Wigner 3j symbols, this is much easier to calculate than 
Gaunt 's formula of the integral. 

The factor rf in (4.4) for j — 2 can be incorporated by the identity (2.14) 
above. 

4.2 Radial Coulomb Integral 

For conciseness of notation, we subsume the "quantum numbers" n£m from (4.1) 
as i, and equivalently for i', i and i' . The radial contribution to the Coulomb 
integral is computationally much more challenging due to the dependence of 
whether £1 < £2 or £1 > £2- Thus, the integrals over £1 and £2 cannot be 
separated any more; instead, we obtain the nested integrals 

Ai(&)M£ira&) Si d ?i d ^ + <«'i ++ h'j) 

(4.6) 
f° r 3,3 G {0: 2} due to the volume element. 

The authors [2] apply an integral transformation (from Ref. [27]) to (4.6) 
and then solve the resulting integral numerically. It consists of an outer integral 
over the product of two functions, which are themselves integrals. Although 
this approach inherently respects the symmetry k -h- k, we haven't found it 
computationally advantageous as compared to solving (4.6) directly, since three 
integrals need to be calculated instead of two. 

In the following, we provide the details of our approach. We employ the 
methods developed in section 3 to evaluate (4.6). As first (and most expensive) 
step, set pa> := (pi +jv)/2 and calculate dw via (3.12) such that 

A, ; (cf)A,, (£) = tf£ (2 Vl x)H^ (2 Pi , x) = H^, {2 Vii , x), x := £ - 1. (4.7) 

Proceed analogously for p- v and <%, . Finally, set p iV ^i '■— {pu> + p^;)/2 and 
calculate coefficients bu>, b^, via (3.13) such that 

H^ i ,(2p ul x) = H^ u/ (2p iilTi ,x) 
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(cquivalently for bjj, ) . In case pu> = pij, , this step can be cut short by simply 
setting bu' := dw and b^, :— <%, . Then, after a change of variables, the 
integral (4.6) (times the normalization factor (r — v)\/(t + v)\ and for j,j = 0) 
equals 

(bu'lBMbrr,), z:=2p uqv (4.8) 

with the matrix B"(z) = ( b v ,Az) ) . defined by 



r.kk 



kk 



Km (*) == j^f f W*)QW + */*) 



#£(z)P^(lH-a;/z)da;dz + /fc-H- A;\. I 4 - 9 ) 



The factors fj and ^ f° r j = 2 or j = 2 in the integral (4.6) can be incorporated 
via equation (2.16), similar to the angular integral. 

Thus, given the matrix B"(z), we have reduced the rather expensive inte- 
gral (4.6) to the simple matrix formula (4.8). To obtain B"(z), we have first 
precomputed the entries (4.9) symbolically in z as detailed below. Still, the 
resulting formulas are quite extensive and preclude a fast numerical evaluation. 
Our remedy consists in a Taylor expansion of (the entries in) B^(z), 

"max / _ \Tl 

B^z)«J2^—^-B^\zo). (4.10) 

We precompute the derivatives B" symbolically (up to n max = 8) and then 
evaluate them at (half)-integers zq = 1, 1.5, 2, 2.5, .... Due to potential numeric 
cancellation effects, we employ high-precision arithmetic for this intermediate 
step. Nevertheless, the entries of the resulting matrices B" {zq) are well- 
behaved and do not increase for higher values of n. These numeric matrices are 
then stored on disk for later usage. 

Error estimation The sampling of half-integer evaluation points ensures that 
for each occurring z, there is a closest z with \z — z \ < 1/4. Thus, a very coarse 
error estimate of the Taylor expansion (4.10) gives an error of 10 -11 , when 
assuming that the individual entries of Br are in the order of 1, independent 
of n. In reality, we observe even better results, up to double floating-point 
precision 10 -16 . 

Until now, we have not yet discussed the truncation error of the Laguerre 
expansions. As illustrated in figure 2 above, we can reach machine precision due 
to the exponential decay. However, the number of required coefficients depends 
on the particular decay parameters. When multiplying two expansions via (3.1), 
these numbers typically add up to give the number of coefficients in the resulting 
expansion. Thus, in our setting, the coefficient vectors ba' and 6^, from the 
formula (4.8) have approximately length 36. 

Cost analysis Summarizing the above steps after precomputation, our al- 
gorithm only needs the numeric matrices Br (z ) from the Taylor expan- 
sion (4.10) as input, instead of the symbolic integrals (4.9). In particular, no 
numeric integration is required. 
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The most expensive remaining step is the Laguerre expansion of the product 
Aj(£)Aj'(£) in (4.7). Assuming that the expansion vectors di and d^ have 
length K and the resulting vector dw length 2K, the operation (3.2) has to be 
performed 2K times, leading to the asymptotic total cost 0(K 3 ). In our setting, 
K is typically equal to 18. Since matrix operations are highly optimized, the 
computation time is in the order of milliseconds on modern PCs. 

Symbolic calculation of the integrals (4.9) In what follows, we reduce (4.9) 
to the integrals (3.16), (3.19) and (3.20) (solved in section 3.3). We focus on 
the relevant cases v = 0, 1, 2, but our approach can easily be extended to higher 
v. In the simplest case r, v — 0, the integrals (3.20) and (4.9) coincide since 
Qo(£) = arcoth(£). For general r, v, our strategy consists of "absorbing" the 
Legendre functions into the Laguerre polynomials from H% and H~ by repeat- 
edly applying equation (2.15) (multiplication by x). 

First, remember that the Legendre function of the second kind splits into 

Qr(0 = ,^l/2 + P r(0 arcoth(0, 

(? - 1) 

where P" (£) is the Legendre function of the first kind and G V T a polynomial of 
order r — 1 + v. 

Consider the case v = 0: Since -P T (£) is actually a polynomial, repeated 
application of (2.15) allows us to write 

k+T 

L k (x)P T {l + x/z) = Y^ a kTt (z)L l (x) (4.11) 

i— k — r 

with some coefficients ak T i(z). Proceeding similarly with G V T (\ + x/z), the inte- 
gral (4.9) becomes as a linear combination of the integrals (3.16) and (3.20). 

The case v = 1 is more involved since -P^(£) is no polynomial any more, but 
consists of the factor (£ 2 — l) 1 / 2 x a polynomial of order r — 1. To circumvent 
this difficulty, we first rewrite 

H l k (x)P^(l + x/z) = —L=H l k {x) ■ ^/27T^P T 1 (l + x/z). 
\J2z + x 

A symbolic Laguerre expansion of the first factor (computed in (3.21), sec- 
tion 3.3) transforms the right hand side to a linear combination of 

H}(x)V2z + xP}(l + x/z). 

Plugging in the definition of H}(x), we obtain 

VWTI) L\ (x) e' x/2 y/x(2z + x) P\ (1 + x/z) . 

poly(a;) 



Since the second part is a polynomial in x, we can proceed as for v = 0. The 
same transformation works for Q V T (Z,) as well, which completes the case v = 1. 
Finally, for v — 2, the decomposition (4.2) times the factor x (from H^(x)) 
reads 



xQl(l+x/z) = G 2 T {l + x/z) + x P 2 {1 + x/ z) &vcot\i{l + x/ z). 
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After a transformation similar to (4.11), we conclude that the integral (4.9) can 
be reduced to a linear combination of the integrals (3.20) and (3.19) (after a 
change of variables x — > x/2 and using the Laguerre product coefficients cf (1) 
in (3.6)). 

Testing the implementation As first check, we calculate the following (homonu- 
clear) Coulomb integral and obtain 

(iplsa^Ua | ^l«r^l*r) = 0.780883 for #=1.4a.u., Z a = Z b = l 

with summation up to r max — 9 in (4.2). This agrees to all digits with the 
tabulated value in Ref. [3]. Similarly, for the heteronuclear case (HeH + molec- 
ular ion with Z a = 2 and Z b = 1), we obtain (i/'iso-^iso- | V'lso-'^iso-) = 1.23207, 
which agrees in 4 digits with the value 1.23225 from Ref. [3]. The discrepancy 
could stem from lower precision arithmetic in [3], or from a different trunca- 
tion of the Neumann expansion. We have verified all digits of our value using 
Mathematica's numeric integration routines to solve the integrals (4.4) and (4.6) 
directly (which is much slower in this case), and by comparing the truncations 
r max = 6, 7, 8, 9 of the Neumann expansion (all agreeing in the first 7 digits). 

For further testing, we have numerically computed the integrals (4.6) with 
several other parameters, and found that the values agree in at least 12 digits. 

5 Application to Diatomic Molecules 

To demonstrate the feasibility of our algorithmic approach, we apply it to the 
diatomic molecules O2 and CO. 

The ./V-electron Hamiltonian for diatomic molecules in atomic units (Born- 
Oppenheimcr approximation) reads 

H = J2(- 1 -A. l -^-^) + J2- + ^, (5.1) 

y\ 2 r *a r lb J f^m R 

where r ia and r^, denote the distances of the i th electron to the fixed nuclei 
at (0,0, =Fi?/2), respectively, Z a ,Zb € N >0 the nuclear charges (as for the sin- 
gle electron Schrodinger equation (2.1)), and nj = \x{ — Xj\ the inter-electron 
distance between electron i and j. The first sum (denoted Hq) contains pre- 
cisely the single-electron Hamiltonian (2.1), the second sum is the inter-electron 
Coulomb repulsion (denoted V ee ), and the last term the repulsion of the nuclei. 
The homonuclear version (Z a = Z b =: Z) describes atomic dimers like hydrogen 
H2 or oxygen O2. 

Analogous to Ref. [12], it is instructive to investigate the limit of large nuclear 
charge Z. (We consider the homonuclear case here for simplicity.) Namely, a 
short calculation shows that if }&(xi,(Ji, ■ ■ ■ ,xn,pn) G Ll((R 3 x{±^}) N ) solves 
the TV-electron Schrodinger equation H^ = E"i> with H defined in (5.1), then 
the rescaled wavefunction 

*(l/!, cn, ...,y N ,<r N ):= Z~ 3N / 2 tffZ-^yx, cti, . . . , Z- l y N , a N ) 

solves 

#0 + ^ + 1)* = -** (5.2) 
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itli 




Ha := , , .^, , ,..._, , , 

3/i - 5/ite 3 



As Z — > oo, 7?o describes two isolated atoms and the electron-electron interac- 
tion \V e e becomes small due to the prefactor -%. Since Hq depends on Z, wc 
cannot repeat the exact same analysis as in Rcf. [12], but our single-electron 
wavefunctions are eigenfunctions of H nevertheless. Thus we expect that our 
calculations match highly-charged (electrically confined) molecular ions well and 
could serve as benchmark for alternative computational approaches. 

To allow for comparison with experimental data, we focus on the paramag- 
netic "triplet" oxygen molecule O2 in the following paragraph, i.e., Z = 8 and 
N = 16. The ground state symmetry is characterized by the molecular symbol 
3 Sg . That is, the spin quantum number equals 1 (hence "triplet"), the angular 
L z momentum quantum number is zero (rotation about internuclear axis), and 
the parity is even. 

The common textbook version of the electronic quantum state reads as fol- 
lows. All molecular orbitals up to ir u (2px,y) are completely filled (see figure 4), 
leaving the two remaining electrons in the "antibonding" tt* (2p xy ) orbitals. 
These two electrons form a spin triplet, hence the paramagnetism. With the 
mapping from figure 4, the electronic configuration corresponds to the following 
Slater determinant: 

*i = |lsa g t| lpcr u n 2scr g n 2p<7 u n ldo-gtt lp(±7r) u U ld(±7r)gt) • 

In what follows, we try to approximate the groundstate energy of O2 via 
the methods from the previous chapters, with summation up to r max = 9 in the 
Neumann expansion (4.2). We include all wavefunctions of the 3 E~ symmetry 
subspace, restricted to configurations with the lscr g ti and lpcr u yi. orbitals com- 
pletely filled, and the occupations of the higher orbitals (up to If a u fi) allowed 
to vary. In our case, this gives 54 wavefunctions, including $1. For example, 
another state in the symmetry subspace reads 

* 2 = |---lp(±7r) u tld(±7r) g n), 

which is the same as ^1 except for half-occupied lp(±7r) u molecular orbitals 
instead of ld(±7r) g . 

Thus, the groundstate energy is the smallest eigenvalue of the 54 x 54 matrix 
(\Pi I H^!j) i , with the Hamiltonian H defined in (5.1). Since the \£^ are exact 
eigenstates of the TV-body Hamiltonian without the inter-electron Coulomb re- 
pulsion V cc , the latter can be regarded as perturbation of (H — V cc ) (see also 
equation (5.2)). 

We use the software toolbox [26, 25] to express (^>i | V 00 ^j) as linear combi- 
nation of Coulomb integral symbols (4.1), after tracing-out the spin variables. 
The symmetry properties (ab | cd) = (cd \ ab) and (ab | cd) = (ba \ dc) simplify 
the resulting expressions. As concrete example, the following off-diagonal ma- 
trix element reads 

(*1 I ^00*2) = (tplp^lpld-Kg I V ; lp(-7r) u V ; ld(-7r)g)-(V'lp7r u '01d(-7r) g I V'lp(-7r) u V'ld 7 r g ) 

Both diagonal entries (^ | V cc ^i) (i = 1, 2) are quite extensive, consisting of 79 
individual Coulomb integrals. 
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Our approach can easily be adapted to other symmetry subspaces. Thus we 
include the experimentally next low- lying symmetry levels 1 A g and ^^ as well 
(see for example Ref. [4] for an overview). 
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Figure 8: Potential energy curves of the O2 molecule (this paper), i.e., low- 
est eigenvalue of the matrix (\l/j | iJ\fj) with the Hamiltonian H from (5.1) 
restricted to the corresponding symmetry subspace. Dots are calculated val- 
ues, and the continuous line a spline interpolation of degree 3. The minimum 
— 133. 689 a. u. of the 3 E~ symmetry level is attained at R — 2.325 a.u. (dotted 
red line). For comparison, the experimental bond length from the literature 
(green line) is also shown. 

The result of our calculations is plotted in figure 8, showing the electronic 
groundstate energy (blue dissociation curve) as well as excited energy levels 
(purple and brown curves) dependent of the nuclear distance R. Our compu- 
tation predicts an optimal bond length R mm = 2.325 a.u. (dotted red line), 
which is quite close to the experimental value from the literature [31, 17], 
i? eX p( 16 02) = 121pm = 2.2866 a.u. (green line). Additionally, we reproduce 
the experimental ordering of the symmetry states. 

Having obtained the groundstate energy, we can calculate the dissociation 
energy O2 — > 2 O by subtracting (2x) the energy of an individual oxygen atom. 
Since the outcome of theoretical calculations depends on the particular model 
(e.g., the Ansatz space of single-electron wavef unctions) , similar models should 
be used for both the O2 molecule and the individual atoms. In our case, a 
close match regarding single atoms is Ref. [12] as already mentioned above. 
Namely, the authors use hydrogen-like wavefunctions as Ansatz space and treat 
the inter-electron Coulomb repulsion as perturbation (similar to the present 
study). Additionally, the electronic configurations match ours in the R — > 00 
limit (available atomic subshells Is, 2s, 2p, with the lowest Is subshell always 
occupied). From [12], -Eo.min = —66.7048 a.u. for the groundstate angular mo- 
mentum/spin symmetry 3 P. Thus, we obtain the dissociation energy 



2Eq, 



Eo 2 



0.278971 a.u. 



(5.3) 



For comparison, the experimental dissociation energy of oxygen (enthalpy change 
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at OK) is -Eo 2 ,ex P diss = 5.1157 eV = 0.1879 a. u. [8], which differs from our cal- 
culated value (5.3) by approximately 50%. The discrepancy is likely due to the 
small dimension of the Ansatz space (number of single-electron wavefunctions, 
up to the 2p subshell in our case). Note that the dissociation energy is 3 or- 
ders of magnitude smaller than the total energy. Thus, subtracting groundstate 
energies requires at least 4 correct decimal digits for just 1 digit of the dissocia- 
tion energy. In any case, our calculated value reproduces the experimental data 
qualitatively correct, in particular the sign (i.e., the fact that O2 binds). 

To provide an example for a heteronuclear molecule, we repeat the analo- 
gous calculations for carbon monoxide CO, i.e., Z a — 8, Z\, = 6 and N = 14. 
Figure 9a shows the resulting ground state dissociation curve with the same 
spheroidal Ansatz space (up to l/crfl) as for oxygen. Notably, the deviation 
between the experimental bond length i? exp ( 12 C 16 0) = 112.8pm (green line, 
[6, 17]) and the calculated minimizer of the curve (dotted red line) is relatively 
large. This is presumably due to the small number of spheroidal basis func- 
tions. Indeed, when including the Sscrti spheroidal orbitals, the minimizer of 
the curve approaches the experimental value (figure 9b). 
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Figure 9: Potential energy curve of the CO molecule restricted to the ground- 
state 1 E + symmetry subspace (this paper), (a) Same basis set as in figure 8; 
(b) additionally including the 3scrt4- spheroidal orbitals. In (b), the minimizer 
of the curve (dotted red line) is closer to the experimental bond length (green 
line). Note that the energy axes are shifted by approximately 3 a.u. 



6 Conclusions and Outlook 

We have developed and implemented an efficient computational framework to 
evaluate the angular and radial Coulomb/exchange integrals in prolate spheroidal 
coordinates by employing Neumann's expansion of 1/ \x — y\ and taking advan- 
tage of symbolic integration as far as possible. The algorithm strongly relies on 
matrix operations to speed up computations. 

A particular advantage of our approach is the universality of the precom- 
puted numeric matrices in (4.10). Once obtained, these matrices can be reused 
for subsequent calculations. 

The application to the oxygen and carbon monoxide molecules shows the fea- 
sibility of our algorithm. We reproduce qualitatively correct energy curves, and 
the calculated bond length and dissociation energy are in reasonable agreement 
with experimental values. 
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A long-term goal of the present paper is a better understanding and quan- 
titative description of atomic interactions and chemical bonds, which could be 
modeled using spheroidal orbitals. To reduce complexity, one could employ the 
well-known hydrogen-like orbitals for the core electrons (close to the nucleus). 
This combination of spheroidal and hydrogen-like orbitals requires proper or- 
thonormalization and the calculation of Coulomb/exchange integrals between 
these different kind of orbitals. Inversing the LCAO Ansatz to approximate the 
spheroidal wavefunctions locally (close to an atomic nucleus) might be feasible 
for these purposes. 

Finally, the algorithm presented here could be combined with established 
computational chemistry methods (like Configuration Interaction or Coupled 
Cluster) in future projects. 
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